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ABSTRACT 

This paper describes a computer program developed for 
structural dynamic analysis of horizontal axis wind 
turbines (HAWTs). It is based on the finite element 
method through its reliance on NASTRAN for the devel- 
opment of mass, stiffness, and damping matrices of the 
tower and rotor, which are treated in NASTRAN as 
separate structures. The tower is modeled in a sta- 
tionary frame and the rotor in one rotating at a 
constant angular velocity. The two structures are 
subsequently joined together (external to NASTRAN) 
using a time-dependent transformation consistent with 
the hub configuration. Aerodynamic loads are computed 
with an established flow model based on strip theory. 
Aeroelastic effects are included by incorporating the 
local velocity and twisting deformation of the blade 
in the load computation. The turbulent nature of the 
wind, both in space and time, is modeled by adding in 
stochastic wind increments. The resulting equations 
of motion are solved in the time domain using the 
implicit Newmark— Beta integrator. Preliminary 
comparisons with data from the Boeing/NASA MOD2 HAWT 
indicate that the code is capable of accurately and 
efficiently predicting the response of HAWTs driven by 
turbulent winds. 

INTRODUCTION 

Throughout the history of the DOE-sponsored horizontal 
axis wind turbine (HAWT) program efforts have been 
undertaken to develop tools for the structural dynamic 
analysis of HAWTs. A number of capabilities have 
emerged, including natural mode and frequency calcu- 
lations with NASTRAN, rigid— rotor aerodynamic load 
codes, dynamic flexible— rotor codes fixed in space at 
the hub, and full dynamic models of the rotor turning 
on the tower . 

The NASTRAN mode and frequency analysis is capable of 
tracking some of the frequencies as they increase with 
increasing rotor speed due to centrifugal stiffening 
but is not adequate for those which are sensitive to 
other rotating coordinate system effects. 

The rigid-rotor codes require the rigid body motion of 
the rotor as input and then compute the aerodynamic 
loads along a blade as it moves through one revolution 
for steady atmospheric conditions. The calculated 
loads are integrated to obtain static section loads 
and moments at any station. Even with this simple 
model, if a reasonable rigid body motion is pre- 
scribed, mean loads are predicted with good accuracy. 
However, the vibratory flapwise loads are generally, 
substantially underpredicted. An example of a rigid- 
rotor aerodynamic load code is the PROP software [Ref. 
l], developed at Oregon State University. 


*Thi s work performed at Sandia National Laboratories 
supported by the U. S. Department of Energy under 
Contract Number DE-AC04-76D00789 . 


The flexible-rotor codes use aerodynamic load models 
which are similiar to those used in the rigid-rotor 
case, but motion of the rotor relative to the hub is 
permitted. Thus the motion of the rotor as well as 
the applied loads are computed. Through the inter- 
dependence of rotor motions and the aerodynamic loads, 
this software accounts for aeroelastic effects. With 
these codes, as before, the mean response of the rotor 
is accurately predicted, but the vibratory response is 
underpredicted. Probably the most widely used of 
these codes is MOSTAB [Ref. 2], which is a derivative 
of a code developed by Paragon Pacific for the dynamic 
analysis of helicopters. 

The full dynamic models add the interactions between 
the tower and the rotor to the f 1 exibl e— rotor software 
described above. Within the confines of small dis- 
placement theory, the rotor is modeled in a rotating 
frame, the tower in a fixed one, and the two struc- 
tures are connected using time-dependent constraints 
or forces. Generally a transient integration tech- 
nique is used to solve the resulting equations of 
motion. Even with the increased level of sophis- 
tication, these codes still underpredict the vibratory 
response. Two examples are the MOSTAS code [Ref. 3], 
which is from the same family as MOSTAB, and DYLOSAT 
[Ref. 4], a proprietary code developed by the Boeing 
Aerospace Company . 

The software described here, named HAWTDYN , is of the 
full dynamic model class. Two features ^ihich set it 
apart from other codes in this class are its use of 
NASTRAN for mass, stiffness and damping matrices, and 
output processing, and the inclusion of stochastic 
wind increments in the aerodynamic load computation. 
Time— dependent constraints, which produces time- 
varying coefficients in the equations of motion, are 
used to connect the rotor to the tower. Aerodynamic 
loads are computed using interference factors pre- 
dicted by the PROP code [Ref. l] for a pre-established 
rotor orientation. Aeroelastic effects are included by 
incorporating the local velocity and twisting defor- 
mation of the blade in the load computation. The 
stochastic wind increments are computed by the method 
outlined in Ref. 5. The resulting equations of motion 
are solved in the time domain using the implicit 
Newmark-Beta integrator. 

Results are presented for a model of the MOD2 wind 
turbine which was designed and fabricated by the 
Boeing Aerospace Company (BAC) . In addition to a fan- 
plot, which shows how the natural frequencies of the 
turbine vary with operating speed, structural load 
time series have been obtained for two stochastic 
winds, one with a mean of 20 mph, and the other 27 
mph . These lime series are reduced and compared with 
field measurements. In order to determine the effect 
of the tower on the structural response, the model was 
modified to fix the rotor hub With this alteration, 
HAWTDYN is consistent with the codes of the flexible- 
rotor class except for its inclusion of stochastic 
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wind effects. Results are also presented for this 
model . 

The following sections contain a description of the 
mathematical model upon which HAWTDYN is based, the 
details of the MODS finite element model, presentation 
and discussion of results, and concluding remarks. 

HAWTDYN MATHEMATICAL MODEL 


If the matrices of Eqn_ (2) are renamed M, C and K, 
and the for ci vector, F, the following equation is 
obtained from the combi nation of Eqns . (2) and (3), 
and premultiplication by A : 

(a t ma ) U + ( A T CA+2A T MA ) u 
+ ( a t ka+a t £\‘+a t ca ) u = a t f. 


Due to its power and versatility in modeling struc- 
tures, the finite element method has been chosen as a 
framework for the derivation of the equations of 
motion for the HAWT . For this derivation, two coor- 
dinate systems are employed in order to represent 
motions throughout the structure as small relative to 
the appropriate frame. Thus the tower is modeled in a 
fixed frame and the rotor in one which rotates at the 
operational speed of the turbine about an axis which 
is fixed in space. The origins of both coordinate 
systems are fixed at the initial hub location. For 
the latter case, rotating frame effects, such as 
Coriolis and centrifugal forces, must be included. 
Considering the tower and rotor as separate struc- 
tures, the equations of motion for each are repre- 
sented below: 

“t «T + C T *T + U T “ F r (1) 

“r ”r + (VM °R + ( V S q) U R “ F c + F s + V 


The transformation matrix, A, only modifies terms in 
the matrices associated with tower or rotor connection 
nodes, and, by judicious selection of the physical 
modeling at these points, certain terms in Eqn . (4) 
can be simplified. For example, if the tower connec- 
tion node possesses only lumped translational mass, 
the terms, A T Ma , AMA , and Am', are rendered inde- 
pendent of time and need only be computed once. More- 
over, if the tower connection node not directly 
involved in any damping, the term A CA also becomes 
time- independent and A CA vanishes. The remaining 
quantity, AKA , will normally be a function of time 
and must be recomputed at each time step. 

Replacing the coefficient matrices of Eqn. (4) by M, C 
and K, the system equations of motion are obtained and 
presented below: 

MU 4 CU + KU = F . (5) 


Here the subscripts T and R refer to the tower and 
rotor respectively The quantities, and S^, which 
derive from rotating coordinate system effects, are 
the Coriolis and softening matrices. The softening 
matrix accounts for changes in the centrifugal force 
that result from the structural deformations. These 
matrices are developed in detail in Ref. 6. On the 
right hand sides of the equations are the applied 
forces, with the subscripts c, g and a referring to 
the centrifugal, gravitational and aerodynamic forces, 
respectively . 

These equations can be combined into one matrix equa- 
tion as f ol lows : 
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Denoting the time-dependent constraint relation which 
connects the rotor to the tower, consistent with the 
hub configuration, as A, the final set of displace- 
ments, velocities and accelerations, U, U and U, can 
be derived from, 


AU, 


AU 4- AU, 


= AU 4- 2AU 4 AU. 


( 3 ) 


Eqn (5) is complicated by the fact that centrifugal 
stiffening, which arises due to the spanwise 
stretching of the blade under tbe action of the cen- 
trifugal loads, must be taken into account. This 
stretching causes the stiffness to be a function of 
the deformation [Bef . 6], necessitating more complex 
solution procedures. To avoid this complexity, the 
stiffness matrix, in Eqn. (l) is modified to be 

coimnensurate with the quasi-static displacement field 
associated with the centrifugal loads. This is accom- 
plished through iterative solution of the following 
equat i on : 

[MMIM ■ MM* N ■ 

The resulting approximate or mean stiffness matrix 
represents the rotor stiffness for the operating speed 
corresponding to the centrifugal loads in Eqn. (6). 
Thus solutions of Eqn. (5) must be interpreted as 
motions about a prestressed state. 

The aerodynamic forces of Eqn. (1) are computed, 
taking into account blade velocities and deformations 
relative to the rotating coordinate system. This 
provides for the representation of aerodynamic stiff- 
ness and damping in the equations of motion. In order 
to compute these forces, a local blade coordinate 
system is defined using instantaneous unit chord and 
span vectors, e and e , which account for initial 
blade position and prefwist, and the local blade 
deformations. The positive senses are from leading to 
trailing edge for chord, and from hub to tip for the 
span. The third instantaneous unit direction is 
identified as the flap vector, e f , and defined by the 
cross product of the chord and span vectors. The 
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relative wind velocity vector, W ,is given by the 
following expression: 


\ = ( 1 - t ) W » + W si"[ 0 R + Ox ( X B W R)i- (?) 


Here, W is the mean axial wind, which can include 
variations due to wind shear and tower shadow; t is 
the velocity reduction factor corresponding to a trim 

solution for the mean wind; W . is the axial stochas- 

si 

tic wind increment computed using the methods des- 
cribed in Kef. 5; D is the operating speed of the 
turbine; and, is the initial local position vector. 
The direction oT lift, e^, is obtained by taking the 
cross product of W and e , and then adjusting the 
sign of the resulting vector so that its dot product 
with e is negative. The direction of drag, e , is 
subsequently computed to be perpendicular to the 
direction of lift and e , this time adjusting the sign 
so that the dot product of the resulting vector with 
e is positive. With these directions defined, the 
angle of attack, and the lift and drag forces per unit 
length are given by, 


c L («) e L ’ ( 0 ) 


In Eqn. (8), p is the air density, a is the length of 
the chord and C and are the coefficients of lift 
and drag respectively. The quantity, W , is the 
magnitude of the component of the relative wind vector 
normal to the span vector, computed as follows: 


tana = 


L = | pa W/ 


D = | pa W n 2 


W = |W xe I (9) 

n | r s j . 

The lift and drag forces are combined with the gravity 
forces to obtain the total force vector per unit 
length. This vector is numerically integrated along 
the length of each blade element, using a Galerkin 
formulation to obtain concentrated nodal forces. 

Time-domain solutions to Eqn. (5) are obtained numeri- 
cally using an implicit integration technique. For 
equations with constant coefficients implicit methods 
are unconditionally stable, which means that the size 
of the time step is only limited by the desired fre- 
quency resolution. The option for larger time steps 
provides a means to analyze structural response to 
stochastic loading, which requires long-time solu- 
tions. The equations of motion for the HAWT contain 
t ime— dependent coefficients, and therefore, uncon- 
ditional stability is not guaranteed. However, cer- 
tain "ad hoc" procedures can be implemented which 
improve the stability and permit reasonably large time 
steps . 


advertised to numerically dampen these modes. This 
may have been caused by the fact that the algorithm is 
not entirely consistent with the Newton method of 
equation solution. Experience has indicated that 
efforts to make the solution procedure conform to this 
method usually produce a stabilizing effect. For 
example, the stability of the response was signifi- 
cantly improved by changing the evaluation of the 
damping term in Eqn. (5) from the beginning to the end 
of the time step. Because of its conformity to the 
Newton method, the Newmark— Beta implicit integration 
scheme [Ref. 8] was the final choice for the solution 
procedure . 

Eqn. (5), discretized in time according to the this 
algorithm, is presented below: 


“"W + ^t+At + W = 

F t+At, 


U t + At =U t + AU \ + At2 |( 

5*) U i + *W), 

(10) 

=*t + At |<^> u t 

+ ^ U t+At] . 


/$ = . 3025 , y = . 6 . 




The final form of the equation is obtained by making 
the substitutions indicated in Eqn. (10) and rear- 
ranging so that only terms associated with ^ 
appear on the left hand side, as follows: 


M + At r C + At 0K t+M U t+At 


—At ( 1— y)C - At 


(H K 


u 

t+At t 


(ID 


+ l ^' AtK t +i t U t 


+ ' K t+At U t + F t+At . 


Even with the provisions described above to stabilize 
the solution procedures, a small amount of growth 
still occurs for some of the HAWT models that have 
been created. Although the origin of this growth, be 
it physical or numerical, has not yet been determined, 
it can be eliminated by the incorporation of struc- 
tural damping (of the order of 5% of critical). 

In order to avoid duplication of such things as devel- 
opment of finite element matrices, solution proce- 
dures, and input and output processing, the MacNeal— 
Schwendler version of NASTRAN was selected as the 
basis for this development. This particular code was 
chosen for several reasons. First, NASTRAN is a 
general purpose finite element code which contains the 
necessary input options required for producing reason- 
ably accurate models of HAWTs . Solution procedures 
are available which provide for inclusion of centri- 
fugal stiffening effects in the rotor. The DMAP 
programning feature of NASTRAN, which allows the user 
to modify the code without dealing with the FORTRAN 
coding, proved to be very helpful even though it was 
not heavily used. The direct matrix input (DMI) 
option, by which matrices can be modified through the 
input data deck was also invaluable. And finally, 
because of NASTRAN’ s widespread use, familiarity with 
its BULK DATA input lends a degree of user-friend- 
liness to the present software. 


The first implicit scheme implemented, the Hilber- 
Hughes algorithm [Ref. 7], exhibited unstable growth 
in the high-frequency response, even though it is 
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NASTRAN 


HAWTDYN 



Figure 1. HAWTDTN dynamic analysis method. 


The relationships between NASTRAN and HAWTDYN are 
displayed in the block diagram of Fig. 1. The tower 
mass, stiffness and damping matrices are developed in 
NASTRAN relative to a fixed coordinate system. The 
rotor is modeled in a rotating frame with the stiff- 
ness aiatrix reflecting the effects of centrifugal 
stiffening. The Coriolis and softening matrices are 
computed external to NASTRAN find included through the 
DMI input option. These two sets of matrices are then 
supplied to HAWTDYN where they are tied together with 
a rotor/tower connection matrix, which models the 
particular hub configuration Aerodynamic loads are 
obtained using the mean windspeed, the stochastic wind 
increments, and the local blade motion. Although an 
active control system has not yet been incorporated 
into the HAWTDYN software, it would also provide an 
input to the aerodynamic load computation. The 
resulting equations of motion are solved using impli- 
cit direct time integration. Computed displacement 
time histories cen be printed, plotted, and, in some 
cases, spectrally analyzed. The NASTRAN code is re- 
entered tor computation oT structural loads and 
stresses . 

The PROP code [Ref. 1] has been incorporated into the 
HAWTDYN software to supply the local free stream wind, 
interference factors associated with a prescribed 
orientation of the rotor . and lift and drag coef- 
ficients This relationship is shown in Fig. 2, along 
with details of the load computation. After the 
relative wind is obtained using the free stream wind 
modified by the interference factor, the stochastic 
wind increments, and the local blade motion, the angle 
of attack is computed and transmitted to PROP for the 
determination of the coefficients of lift and drag. 



Figure 2. Aerodynamic load computation in HAWTDYN. 


DESCRIPTION OF THE M0D2 FINITE ELEMENT MODEL 

For an initial test of the performance of HAWTDYN a 
finite element model of the M0D2 was created. The 
M0D2 was chosen because of the general interest in 
that turbine and the availability of experimental 
data. The model was actually developed using a 
NASTRAfT BULK DATA deck assembled at the Boeing 
Aerospace Company (BAC) in June of 1981. The BAC 
model consists of 15 nodes per blade and 14 tower 
nodes. For HAWTDYN this model was reduced using the 
NASTRAN ASET option to 5 tower nodes and 5 nodes per 
blade. This reduced model is shown in Fig. 3. Blade 
stations 370 and 1164 are indicated in this figure for 
future reference. 
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Figure 3. MOD2 finite element model. 


To examine the adequacy of this model, predicted 
natural frequencies for the rotor parked in the ver- 
tical position are compared to measured results taken 
from Ref. 9- The table below shows frequencies pre- 
dicted using the original BAC model and the reduced 
HAWTDYN model, and experimental results for several of 
the lower frequency modes. These frequencies are 
normalized by the operating speed of the rotor, which 
is 17.5 rpm. The column labeled percentage error 
corresponds to the HAWTDYN model predictions relative 
to data. 


The tip pitch control is not modeled in an active 
sense in HAWTDYN , but rather the tip pitch is preset 
for use in the aerodynamic load computation. For 
structural purposes the nominal pitch configuration is 
used in all cases. The drive train is modeled with a 
spring and damper attached between the hub and 
nacelle. The damping in the actual hardware, provided 
by the tip control, is not included in HAWTDYN, but 
approximated by setting the drive train damper at 16% 
of critical. The model includes a yaw spring and 
lateral tower damping. The values for the tower 
damping are set at 4% of critical for side-to-side 
motion, and 1% for fore-aft, consistent with measured 
resul ts [Ref . 9] . 

The rotor/tower connection matrix for the M0D2 is 
shown in Eqn. (12). In this equation, sn represents 
sin(Ot) and cs , cos(Ot). The upper case XYZ subscripts 
correspond to the fixed coordinate system, the 
orientation of which is shown in Fig. 3, and the lower 
case ones, the rotating frame. Initially, for the 
blade in the vertical position, these two systems 
coincide, with the origins of each fixed at the hub. 
The rotating system turns in a positive sense about 
the Z axis. The U's and 8’s represent displacements 
along and rotations about the respective axes. 


Generally, when experimental data is available, rea- 
sonably accurate models can be created by modifying 
system parameters in a logical fashion such that the 
errors in the frequency predictions are 5% or less. 
Most of the errors in Table 1 are of this order. 
However, for two of the modes, the symmetric flap and 
the symmetric chord, the errors are particularly high, 
especially considering their importance in the struc- 
tural response of the rotor. Although not pursued 
here, some effort at fine timing is definitely indi- 
cated for a more accurate M0D2 model . 

A fanplot for the M0D2 has been constructed by taking 
fast Fourier transforms on predicted displacement 
histories associated with the f ree— vibrat i ons of the 
turbine. The natural frequencies in this fanplot, 
which is presented in Fig. 4, are normalized to the 
operating speed of the rotor (17.5 rpm). This figure 
demonstrates how the frequencies vary with rpm, and 
indicates their values at the operating speed. The 
frequency corresponding to the synraetric flapwise 
bending mode is of special significance because of its 
proximity to 4/rev. The nearness of this frequency to 
that excitation indicates a possibility for larger- 
than— expected response due to dynamic amplification. 
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Table 1. Comparison of Predicted and Measured 
Frequencies for the Parked M0D2 Rotor 


Mode Freq (/rev) 

Meas . 

Freq (/rev) 
BAC Model 

Freq (/rev) 
HAWTDYN 

% Error 

Drive 

Train 

.45 

.42 

,42 

6.6 

Tower 
Fore/Af t 

1.23 

1.23 

1.24 

.8 

Tower 

Lateral 

1.28 

1.25 

1.25 

2.3 

Symme trie 
Flap 

3.30 

3.77 

3.76 

13.9 

Synme trie 
Chord 

6. 17 

6.79 

6.69 

8.4 

Ant i sym 
Flap 

6.55 

7.03 

7.03 

7.3 

Nacel 1 e 
Pitch 

8.23 

8.71 

8.74 

6.2 

2nd sym 
Flap 

9.60 

10,01 

10.01 

4.3 



Figure 4 Predicted fanplot for the M0D2 . 


PRESENTATION AND DISCUSSION OF RESULTS 

In this section to some qualitative observations are 
discussed, and the predicted response of the M0D2 is 
presented for windspeeds of 20 and 27 mph . In both 
cases a turbulent wind increment history is included. 
For purposes of comparison, results are shown for 
steady winds at the same windspeeds. Computations are 
also presented for a configuration of the MOD2 where 
the hub is fixed in space. 

Two qualitative aspects of the results are mentioned 
here to promote confidence in the HAWTDYN software. 
First, the aerodynamic damping in the fore-aft direc- 
tion is definitely present in the solutions and seems 
to be on the order of 5% of critical. And second, 


HAWTDYN seems to be correctly predicting the response 
of the rotor to wind shear in that computed displace- 
ments indicate a slight turning of the rotor out of 
the wind about a vertical axis. This motion produces 
a more uniform relative velocity vector with respect 
to the angular position of the rotor, and tends to 
neutralize the effect of wind shear. These results 
are qualitatively consistent with observed behavior. 

For the forced response of the turbine, gravity and 
wind loading are applied. The wind loads correspond 
to a wind shear resulting from a surface roughness 
factor of .25, a value consistent with the M0D2 site 
at Goodnoe Hills. In addition, the turbulence inten- 
sity is set at 20%, also representative of the site. 
Tower shadow for a width of 15 feet is included 
through an option in the PROP code. The steady compo- 
nent of the relative wind velocity vector is computed 
at each Gaussian integration point along the blade 
(two integration points per element). The stochastic 
increments are calculated at three stations per blade 
(stations 874, 1329, and 1655) and linearly inter- 
polated or extrapolated to the Gauss points. Fig. 5 
shows the stochastic wind increments for all six 
stations on the rotor, corresponding to the mean 
windspeed of 27 mph. At any particular time the width 
of this band of curves is an indication of the varia- 
bility of the wind across the rotor. 



Tifc rsco 


Figure 5. Stochastic wind increments for a 27 mph 
mean wind at the six M0D2 rotor stations, 
three per blade . 


The predicted edgewise response to this wind at sta- 
tion 370 is presented in Fig. 6. For the first three 
revolutions on the curve, only gravity and steady wind 
loads are acting. Thus these three cycles represent 
the steady response of the turbine. Over the forth 
cycle the stochastic wind increments are gradually 
included. The total response of the rotor to all 
loadings is shown from revolution 5 through 23. For 
all sections of the curve the response is predomin- 
ately 1/rev with only slight variations in amplitude. 
This indicates that the edgewise moments are dominated 
by the gravity loading, and are not significantly 
affected by the turbulence in the wind. 
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TIME (ROTOR REU0LUTIDN5) 


Figure 6. Predicted chordwise bending moment for a 
turbulent wind with a 27 mph mean at 
station 370 of the M0D2 blade. 



TIME (ROTOR REUOLUTIONS ) 


Figure 8. Predicted flapwise bending moment filtered 
using a highpass filter with a cut-off 
frequency of .25/rev. 


The flapwise response at station 370 is shown in 
Fig. 7. Again the first three cycles represent the 
steady response and cycles 5 through 23 display the 
total response. In this case the differences between 
the steady and total responses are quite dramatic. In 
addition to a roughly four-fold increase in the cyclic 
amplitude, the frequency content changes from 2/rev to 
predominately 4/rev. Thus, for the flapwise moments, 
the stochastic wind loading dominates the response, in 
contrast to the edgewise case. 



Figure 7. Predicted flapwise bending moment for a 
turbulent wind with a 27 mph mean at 
station 370 of the MOD2 blade. 


Two additional forms of data reduction are done on the 
curve of Fig. 8. First it is transformed to determine 
its frequency content using the fast Fourier trans- 
form. The transformed curve, which is shown in Fig. 

9, clearly indicates the dominance of the 4/rev compo- 
nent of the response. This dominance is the result of 
the nearness of the frequency of the symmetric flap- 
wise bending mode to 4/rev at the operating speed. 
Consistent with the steady wind response, a signi- 
ficant 2/rev component is also present in the total 
response . 
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Figure 9. Fast Fourier transform of the filtered 
flapwise bending moment predictions. 


In order to compare with experimental data, the same 
types of data reduction that are done on the field 
measurements must be done on the predicted results. To 
this end, the curve of Fig. 7 is first truncated to 
delete the first four cycles, f i 1 tered wi th highpass 
filter to eliminate frequencies below .25/rev and 
finally truncated again to delete spurious results 
near the end of the record caused by end— effect pro- 
blems associated with the filtering. The end product 
is presented in Fig. 8. 


The second form of data reduction consists of cycle 
counting to obtain a 50 percentile cyclic value. The 
first step of this procedure involves tabulating all 
the maximum and minimum values between the zero 
crossings of Fig. 8. The 50 percentile cyclic value 
is then obtained by averaging the absolute values of 
each of the tabulated results. 
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In Fig. 10(b) the 50 percentile cyclic flapwise moment 
at station 370 is displayed as a function of wind- 
speed. The experimental data denoted in this figure 
by the "plus' 1 signs was collected from the M0D2 
cluster at Goodnoe Hills in July of 1983. Predictions 
for two mean windspeeds „ 20 and 27 mph, are indicated 
on the plot by the solid circles. In both cases the 
predictions fall within the scatter band of the data. 
The solid triangle corresponds to the predicted value 
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Figure 10. 
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Comparison of predictions, denoted by the 
solid circles and triangle, to data, 
denoted by the plus signs, at station 370 
of the MOD2 blade, (a) mean flapwise bend- 
ing moment, (b) cyclic flapwise bending 
moment . 


for a steady 27 mph wind. The importance of including 
turbulent wind effects in the structural dynamic 
analysis of HAWTs is clearly demonstrated by this 
figure. 

Fig. 10(a) shows the corresponding mean value of the 
flapwise moment at station 370. As in the cyclic 
case, the predictions fall within the scatter band of 
the data. Although not displayed on the plot, the 
steady wind prediction at 27 mph coincides with the 
value shown for the total response. This indicates 


that the inclusion of stochastic wind effects does not 
seem to be critical for accurately computing mean 
flapwise moments. 

For station 1164, results similar to those of Fig- 10 
are shown in Fig. 11. Again the predictions are 
generally within the data scatter band, with the mean 
flapwise moment at the high rim of the band This 
slight overprediction is probably caused by the fact 
that the continuous blade loads are integrated to form 
concentrated nodal forces. In the static case for a 
uniform load, this discretization of the load produces 
moments which are correct at the node points, but 
overpredicted everywhere in between. 
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Figure 11. Comparison of predictions, denoted by the 
solid circles and triangle, to data, 
denoted by the plus signs, at station 1164 
of the MOD 2 blade, (a) mean flapwise bend- 
ing moment, (b) cyclic flapwise bending 
moment . 


To examine the role that the tower plays in the 
response of the M0D2 , the hub of the rotor was 
constrained so that it could not translate. This did 
not compromise its ability to teeter however. For the 
same turbulent wind with the 27 mph mean, results show 
a reduction In the 50 percentile cyclic flapwise 
moment of approximately 10 percent. Thus, for the 
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M0D2, it may be possible to eliminate the tower from 
the analysis without significantly degrading the 
accuracy of the results. This elimination 
considerably simplifies the analysis procedures. 
However, this observation for the M0D2 cannot be 
generalized to all HAWTs. The critical issue in 
excluding the tower is the degree to which its pres- 
ence modifies the natural frequencies of the rotor. 

The computer resources required by HAWTDYN are modest. 
The M0D2 model , which contains 67 degrees of freedom, 
requires 240 cp seconds on the CRAY computer to obtain 
solutions out to 80 seconds. The relatively large 
time step of .008 seconds, made possible by the use of 
an implicit integration method, was used for these 
calculations. Although tests were conducted to estab- 
lish the accuracy of this time step, no attempts were 
made to find the largest possible time step consistent 
with accuracy and stability restrictions. Thus, it 
may be possible to reduce the cp time below the value 
reported here. 

CONCLUSIONS 

In the design and analysis of dynamic systems three 
areas of concern are routinely addressed: the natural 

frequencies of the system, the excitation frequencies, 
and the ability of the distributed forces to excite 
the natural modes. For HAWTs, the task of addressing 
these areas of concern is not routine. The identi- 
fication of the natural frequencies of the system is 
complicated by the fact that the frequencies must be 
obtained for the turbine in its operational config- 
uration. From the analysis point of view, this en- 
tails connecting the rotor, which moves relative to a 
rotating frame, to the tower which moves relative to a 
stationary one. This precludes the use of standard 
solution techniques for obtaining the natural frequen- 
cies and modes of the structure. The determination of 
the excitation frequencies and the spatial distri- 
bution of the forces is more involved because of the 
turbulent nature of the wind. Not only is the turbu- 
lence difficult to model, but sophisticated methods 
for predicting the resultant loads have not been 
developed. However, with proper attention to creating 
structural and aerodynamic load models that contain 
the major physical aspects of the problem, reasonably 
accurate results can be obtained. 

In the present case, the HAWTDYN software has produced 
accurate predictions for the M0D2 turbine, for two 
windspeeds. These results are of a preliminary nature 
and should be viewed as such until more of the 
validation process has been completed. This process 
includes making predictions for several existing 
turbines at a number of wind conditions, and comparing 
the results with available experimental data. In 
order to generate more confidence in the present M0D2 
finite element model, it should be upgraded so that 
the computed natural frequencies are in better agree- 
ment with measured ones and some development effort 
should be expended in HAWTDYN to properly model the 
tip control mechanism. 

Even though some aspects of the modeling are crude, 
HAWTDYN has produced some promising results. The 
inclusion of turbulent wind effects has dramatically 
influenced the MOD2 response predictions, bringing 
them into agreement with measured data. The use of 
implicit integration methods has posed no insur- 
mountable problems and has made it possible to obtain 
the long-time solutions required to analyze HAWTs 
driven by turbulent winds. After the validation pro- 


cess has been completed, it is expected that HAWTDYN 
will be suitable for accurately evaluating the struc- 
tural response of alternate HAWT designs. 
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